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I.  REPETITA  JUVANT 
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One-dimensional  flows,  in  the  present  paper,  are  unsteady  flows 
which  depend  on  one  space  variable  only.  Therefore,  unsteady  flows 
whose  properties  are  uniform  along  planes,  cylinders  or  spheres,  as  well 
as  flow  in  ducts  of  variable  cross-section,  treated  within  the  framework 
of  the  well-known  quasi-one-dimensional  approximation,  all  belong  to  that 
category.  Despite  their  crudity,  one-dimensional  concepts  find  useful 
applications  in  a  large  number  of  problems,  such  as  shock-tube  analysis, 
intake  and  exhaust  unsteady  behavior,  preliminary  wind-tunnel  design, 
explosions,  stellar  evolution,  vehicle-in-fcube  performances,  etc.  In 
addition,  one-dimensional  concepts  can  be  used,  by  analogy,  in  steady, 
supersonic  flows  depending  on  two  space  variables  (that  is,  two-dimen¬ 
sional  or  axially  symmetric) .  From  the  viewpoint  of  research,  in  gas 
dynamics  as  well  as  in  numerical  analysis,  one-dimensional  problems  pro¬ 
vide  basic  examples  in  the  simplest  possible  way,  clear-cut  ideas  and 
accuracy  tests. 

Has  a  problem  with  so  many  positive  features  and  so  rich  in  practical 
applications  been  formulated  in  a  general,  exhaustive  way?  Is  a  computer 
program  available,  easy  to  use,  general,  safe,  accurate  and  fast? 

Judging  by  the  existing  literature  and  the  requests  from  the  industry, 
the  answer  seems  to  be  on  the  negative. 

The  reason  is  that  one-dimensional  problems  are  not  simple  at  all. 

As  a  consequence  of  the  nonlinearity  of  the  governing  equations,  discon¬ 
tinuities  in  the  physical  parameters  and  their  derivatives  appear,  most 
disturbingly  in  the  form  of  shock  waves  and  interfaces  (contact  discon¬ 
tinuities)  .  Such  discontinuities  interact  with  each  other,  generating 
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more  discontinuities,  and  reflect  on  boundaries.  In  general,  flows  which 
start  with  a  continuous  distribution  of  parameters,  quickly  develop  dis¬ 
continuities.  A  cross-section  of  the  flow  at  a  constant  time  appears  as 
a  complicated  distribution  of  discontinuities,  separating  diminutive 
regions  of  continuous  flow. 

To  make  the  point  bluntly:  Rather  than  as  a  continuous  distribution 
of  parameters,  occasionally  interrupted  by  discontinuities,  a  one-dimen¬ 
sional  flow  can  be  conceived  as  a  pattern  of  discontinuities,  separating 
regions  of  almost  uniform  flow.  Bearing  this  in  mind,  a  computer  program 
should  emphasize  the  handling  of  discontinuities  and  their  interactions, 
rather  than  the  calculation  of  continuous  regions. 

The  trend  in  the  last  two  decades  has  been  in  the  opposite  direction. 

It  is,  by  now,  so  deeply  rooted  that  few  words  to  explain  how  the  shift 

* 

began  may  hopefully  help  to  redirect  our  efforts  onto  what  I  believe  is 
the  right  track.  The  first  attempts  to  calculate  one-dimensional  flows 
were  based  on  the  method  of  characteristics,  a  convenient  technique  for 
hand  computations  and  graphical  analysis  as  long  as  the  flows  are  contin¬ 
uous,  homoentropic  and  strictly  one-dimensional.  In  the  presence  of  dis¬ 
continuities,  the  method  of  characteristics  has  to  be  modified.  The 
latter  can  be  used  in  each  region  of  continuous  flow,  with  the  additional 
complication  of  a  generally  variable  entropy;  the  discontinuities  must  be 
treated  by  different  techniques.  Hand  computations  become  prohibitively 
cumbersome;  interpreting  the  method  in  a  computer  program  is  a  formidable 
exercise  in  logic.  Not  surprisingly,  scientists  interested  in  fluid 

* 

mechanics  shied  away  from  such  problems,  whose  solution  seldom  contributes 
to  deepening  one's  physical  knowledge.  On  the  other  hand,  skilled 
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programmers  who  may  enjoy  a  challenge  in  logic  cannot  undertake  the  task 
without  a  gas  dynamicist ' s  guidance.  In  1950,  von  Neumann  and  Richtmyer 
came  along  with  a  brilliant  idea\  which  seemed  to  sweep  out  all  logical 
difficulties  by  a  single  stroke,  making  gas  dynamical  computations  an 
effortless  routine,  to  be  stolidly  performed  by  the  high-speed  computer, 
over  and  over  again,  invariably,  at  all  the  points  to  be  evaluated.  In 
the  scientific  community,  such  a  technique  soon  became  known  as  "brute 
force",  the  derogatory  tinge  being  intentional.  The  industry,  however, 
welcomed  it  as  a  simple  way  of  obtaining  results,  with  a  minimal  contri¬ 
bution  of  specialized  manpower,  and  of  replacing  expensive  and  delicate 
experiments  by  inexpensive  and  safe  computations.  The  optimistic  term, 
"electronic  wind  tunnel",  was  coined  for  the  high-speed  computer  working 
on  a  gas  dynamical  problem. 

It  turned  out  that  the  computations  were  neither  safe  nor  inexpensive. 
One-dimensional  problems  clearly  show  why  the  method  fails.  To  eliminate 
discontinuities,  one  has  to  smear  them  out  by  introducing  an  artificial 
viscosity  into  the  finite-difference  equations.  If  the  number  of  mesh 
points  is  too  small,  the  artificial  viscosity  is  too  high.  Thus,  dis¬ 
continuities  diffuse  over  too  wide  a  region.  When  more  than  one  discon¬ 
tinuity  exist,  their  broadened  counterparts  tend  to  overlap  and  the  inter¬ 
mediate  region  is  completely  defaced.  A  possible  remedy  consists  of 
increasing  the  number  of  mesh  points.  The  points  must  be  evenly  distrib¬ 
uted  all  over  the  region  to  be  computed  since  (according  to  the  "brute 
force"  approach)  one  does  not  know  where  the  discontinuities  are  going  to 
be  at  any  given  time;  and  the  distance  between  mesh  points  has  to  be  made 
really  small.  Computations  which  may  claim  to  be  safe  (that  is,  accurate) 
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are  catastrophically  expensive2* 


3 


FIG.  1.  COMPARISON  BETWEEN  SOLUTIONS  OBTAINED  BY 

CHARACTERISTIC  AND  MESH  CODES  (FROM  REF.  4) 


A  couple  of  examples  will  help  make  the  point  clear: 

1)  See  Fig.  6  of  Ref.  4  (reproduced  here  as  Fig.  1) ,  where  the 

/ 

dotted  line  plots  the  results  of  a  brute  force  method.  To  quote  the 

authors  of  Ref.  4:  "This  calculation  used  470  points  and  took  three  times 

* 

as  long  as  the  characteristic  calculation  (15  minutes  instead  of  5  ) .  It 
will  be  seen  that  there  is  a  larger  uncertainty  in  the  shock  positions 
and  strengths.  Smoother  profiles  can  be  obtained  by  using  different 
artificial  viscosities,  but  at  the  expense  of  further  shock  broadening  and 

* 

The  computer  used  in  the  test  is  not  mentioned. 
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a  corresponding  increase  in  uncertainty."  I  agree?  I  would  like  to 
anticipate,  though,  that  by  using  the  method  explained  in  the  present 
paper  a  maximum  of  35  points  could  give  accurate  results,  without  certain 
complications  inherent  in  the  method  of  characteristics.  The  reduction  in 
computational  time  would  thus  be  by  a  factor  of  one  hundred  or  more. 

2)  Ref.  5,  instead,  is  a  typical  application  of  brute  force  methods 
with  no  imaginative  criticism.  After  stating  that  "the  availability  of 
larger,  faster  computing  machines  now  allows  the  complete  governing  non¬ 
linear  equations  to  be  solved  by  conceptually  simple  and  more  generally 
applicable  explicit  finite  difference  methods,"  the  author  proceeds  to 
apply  the  obsolete  Lax  scheme,  notorious  for  its  low  degree  of  accuracy. 

The  ratio  of  (CDC  6600)  computer  time  to  real  time  is,  according  to  Ref.  5, 
about  6000  for  relatively  simple  problems  (constant  area  duct,  a  single 
shock) ,  and  should  be  multiplied  by  4  for  more  complicated  problems.  By 
the  method  of  the  present  paper,  the  example  which  requires  200  mesh  points 
in  Ref.  5  could  be  analyzed  with  about  20  points,  thus  reducing  the  compu¬ 
tational  time  by  a  factor  of  100.  The  ratio  of  computer  time  to  real 
time  would  then  be  only  60,  and  the  accuracy  would  be  improved  as  well. 

Note  that  such  a  ratio  is  still  unacceptable  for  many  practical  applica¬ 
tions.  However,  a  machine  ten  times  as  faster  as  the  CDC  6600  would  be 
sufficient  to  make  the  computational  time  acceptable,  and  this  is  within 
our  immediate  reach.  It  should  also  be  noted  that,  for  problems  compli¬ 
cated  by  the  presence  of  many  discontinuities,  the  total  number  of  points 
remains  substantially  unchanged,  if  the  method  of  the  present  paper  is 
used,  whereas  it  is  multiplied  by  a  sizeable  factor  (to  be  squared  to  eval¬ 
uate  the  increase  in  computational  time)  if  the  method  of  Ref.  5  is  used. 
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Along  the  same  line  of  thought,  we  find  in  the  more  recent  litera¬ 
ture  methods  which  seem  capable  of  a  greater  accuracy  than  Lax's  scheme, 
claiming  not  to  use  artificial  viscosity  but,  nevertheless,  to  be  able  to 
find  and  handle  discontinuities  without  explicit  provisions.  One  of  such 
methods  has  been  authoritatively  supported  in  connection  with  other  prob¬ 
lems  of  a  similar  nature^.  The  technique  has  been  labeled  "shock 
capturing"  for  its  purported  ability  to  build  up  sharp  transitions, 
roughly  equivalent  to  shocks,  on  about  three  mesh  points.  In  problems  as 
complicated  as  the  ones  discussed  in  the  present  paper  where,  as  we  will 
see,  regions  described  by  only  three  points  between  discontinuities  exist, 
one  may  wonder  how  the  two  limiting  discontinuities  would  appear. 

We  are  back  again  to  a  need  for  more  points!  However,  regardless  of 
the  number  of  points,  the  method  seems  to  be  inconsistent  because,  if 
viscosity  is  completely  eliminated,  the  equations  of  motion  are  not  prop¬ 
erly  used  in  the  transition  substituting  for  the  discontinuity  (see  Ref.  2, 
pages  19  and  54).  In  the  same  Ref.  2,  page  50,  I  have  mentioned  McCormack's 
scheme  as  a  fairly  good  second  order  method,  and  elsewhere  I  have  repeatedly 
pointed  out  some  of  its  practical  coding  advantages  in  complicated,  multi¬ 
dimensional  problems.  I  cannot  share,  though,  the  optimism  of  the  authors 
of  Ref.  6  in  saying  that  McCormack's  scheme  is  a  "shock  capturing"  one, 
more  than  other  schemes  and  within  acceptable  time  limitations. 

To  submit  my  theoretical  arguments  of  Ref.  2,  page  54  and  Lomax  and 
Kutler's  opposite  statements  to  a  practical  test,  we  may  resume  the  cal¬ 
culation  which  in  Ref.  2  was  halted  soon  before  the  appearance  of  an 
imbedded  shock.  According  to  the  theoretical  analysis  of  Ref.  7,  page  38, 
the  shock  should  form  at  t=.844,  x=.8415.  Results  obtained  by  using 
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McCormack's  scheme,  with  the  original  equations  either  in  conservation 
form  or  not,  and  Ax  equal  to  .05,  .025,  .0125,  .00625  show  no  differences 
from  each  other  as  long  as  t  does  not  exceed  .8. 

If  t  increases  beyond  .8,  the  results  obtained  when  McCormack's 
scheme  is  applied  to  the  non-conservative  form  of  the  equations  oscillate 
in  the  region  where  the  shock  should  appear.  The  oscillating  pattern 
worsens  by  refining  the  mesh,  a  feature  which  makes  any  attempt  to  better 
accuracy  hopeless. 

Better  results,  on  the  high  pressure  side  of  the  shock,  are  obtained 
if  the  equations  are  recast  in  conservation  form.  The  low  pressure  side, 
however,  does  not  improve.  Fig.  2  shows  the  velocity  distribution  at 
t=.901,  as  computed  by  the  McCormack  scheme  with  the  equations  in  con¬ 
servation  form,  with  Ax=. 001667  (360  points),  and  with  Ax=. 013333  (45 
points) ,  together  with  the  velocity  distribution  at  the  same  time,  com¬ 
puted  by  the  method  of  the  present  paper,  using  a  total  of  20  points  (Ax 
equals  .03  before  the  formation  of  the  shock  and  it  maintains  a  value  of 
the  same  order  of  magnitude  after  the  shock  is  formed) .  The  shock  loca¬ 
tion,  properly  found  by  the  present  technique  with  20  points,  is  detected 
by  the  "shock  capturing"  scheme  if  360  points  are  used.  The  low  pressure 
side  distribution  is  misrepresented  even  with  360  points  (note  that  in 
this  problem  the  velocity  in  front  of  the  shock  is  small  but  positive,  as 
the  exact  solution  and  our  results  show) .  One  can  anticipate  some  embar¬ 
rassment,  should  another  shock,  proceeding  from  right  to  left,  impinge  on 
the  former. 

Anyway,  the  ratio  of  computational  time  between  the  360  point  case 
and  the  20  point  case  is  329  (~-  x  *  the  second  fraction  being  the 
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2. 

- McCORMACK  SCHEME,  360  POINTS 

McCORMACK  SCHEME,  45  POINTS 
— O-  PRESENT  TECHNIQUE,  20  POINTS 


FIG.  2.  EXPLICIT  SHOCK  VS.  SHOCK-CAPTURING  RESULTS 
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ratio  of  time  steps  to  reach  t=.901) .  With  such  a  ratio,  assuming  that 
the  present  technique  can  handle  a  problem  at  a  ratio  of  computational 
time  to  real  time  equal  to  40,  the  same  ratio  for  the  shock- capturing 
technique  becomes  13160,  twice  as  larger  as  the  one  mentioned  in  Fef.  5. 

II.  A  MARCH-ON  TECHNIQUE  EMPHASIZING  THE  ROLE  OF  DISCONTINUITIES 

The  present  paper  introduces  a  computational  technique  wnJ  Lakes 
full  advantage  of  explicit  treatment  of  discontinuities.  Although  the 
current  version  of  the  code  uses  the  concept  of  characteristics  occa¬ 
sionally,  its  general  logic  is  not  based  on  the  method  of  characteristics, 
as  in  Ref.  4.  The  belief  that  only  if  the  method  of  characteristics  is 
used,  discontinuities  can  be  treated  explicitly  is  unjustified  and  only 
attributable  to  the  historical,  parallel  development  of  the  method  of 
characteristics  and  finite-difference  schemes  as  outlined  above.  The 
physical  problem  is  sufficiently  complicated  per  se,  not  to  ask  for  the 
unnecessary  logical  complications  brought  in  by  the  method  of  character¬ 
istics.  Moreover,  it  must  be  noted  that  the  latter  does  not  offer  any 
advantage  as  far  as  accuracy  and  computational  speed  are  concerned. 

The  additional  complications  mentioned  above  proceed  from  the  in¬ 
trinsic  irregularity  of  a  characteristic  network.  In  addition  to  making 
plotting  very  cumbersome,  the  characteristic  pattern  can  be  so  different 
in  different  regions  of  the  (x,t)  plane  that,  in  the  course  of  computing, 
uninteresting  regions  or  regions  where  the  computation  may  fail  for  lack 
of  accuracy  are  reached  when  the  region  of  interest  has  not  been  fully 

explored  yat.  Finally,  if  the  computation  proceeds  building  up  a  right¬ 
running  characteristic,  say,  a  left-running  shock  can  be  evaluated  without 
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excessive  additional  labor;  a  right-running  shock, however,  cannot  be 
evaluated  without  first  determining  a  sufficiently  wide,  shockless  region 
across  which  the  shock  has  to  be  laid,  thus  canceling  the  portion  of  the 
computed  region  on  its  high  pressure  side.  Logical  difficulties  arise 
and  grow  bigger  and  bigger  as  tne  discontinuities  interact  with  each 
other.  In  addition,  time  is  lost  in  computing  parts  which  are  succes¬ 
sively  erased.  Finally,  since  the  latter  have  no  physical  meaning,  out- 
of-range  numberscan  be  produced  and  the  computation  halted. 

None  of  the  above  difficulties  exists  if  one  proceeds  from  a  given 
time  to  a  successive  time,  advancing  by  the  same  amount  At  at  all  points. 
This  is  easily  and  accurately  done  by  using  a  suitable  finite-difference 
scheme  at  all  points  where  the  physical  parameters  are  continuous  and 
differentiable,  and  their  first  derivatives  are  continuous. 

However,  according  to  the  spirit  of  the  present  technique,  I  will 
postpone  the  discussion  of  the  finite-difference  scheme  and  discuss  first 
the  physical  features  of  discontinuities  and  their  interaction. 

III.  SHOCKS,  CONTACT  DISCONTINUITIES  AND  GRADIENT  DISCONTINUITIES 

Consider  the  flow  described  by  three  physical  parameters,  the 
particle  velocity,  u,  the  logarithm  of  pressure,  P,  and  the  entropy,  S. 

A  shock  is  a  locus  of  discontinuities  for  u,  P  and  S.  It  travels  at 
a  speed,  W,  which,  relatively  to  the  moving  particles,  is  always  greater 
than  the  speed  of  sound  on  the  low  pressure  side.  Seven  equations  are 
needed  to  determine  the  values  of  u,  P  and  S,  on  either  side  of  the  shock, 


A  contact  discontinuity  is  a  locus  of  discontinuities  for  S  alone. 

It  travels  with  the  particles,  that  is,  at  a  speed  equal  to  u.  Therefore, 
four  equations  are  needed  to  determine  u,  P,  and  the  two  values  of  S  at 
either  side  of  the  discontinuity. 

A  gradient  discontinuity  is  a  locus  of  discontinuities  for  the  first 
derivatives  of  u,  P  and  S.  It  travels  as  a  characteristic.  At  a  point 
located  on  a  gradient  discontinuity  there  are  no  more  unknowns  than  at  an 
ordinary  point.  However,  one  cannot  replace  a  derivative  by  a  finite 
difference  involving  points  at  both  sides  of  the  line. 

If  each  discontinuity  is  considered  as  a  boundary  point  between  two 
regions  of  continuous  flow,  we  may  count  such  a  point  twice,  once  as 
being  the  last  point  in  the  region  at  the  left  and  again  as  being  the 
first  point  in  the  region  at  the  right. 

Suppose  the  discontinuity  is  a  shock.  Since  the  flow  in  the  low 
pressure  region,  relative  to  the  shock,  is  supersonic,  the  flow  values 
at  time  t  are  necessary  and  sufficient  to  determine  the  shock  point  on 
the  same  side  at  t+At,  provided  that  its  location  is  known.  Both  (u+a) 
and  (u-a)  characteristics,  indeed,  reach  the  shock  from  the  low  pressure 
side,  and  so  does  the  particle  path.  Any  three  independent  equations 
(such  as  the  two  compatibility  equations  along  characteristics  and  the 
equation  expressing  the  constancy  of  S  along  a  particle  path,  or  Euler's 
equations  in  their  original  form)  provide  the  values  of  P,  u,  and  S  at 
the  shock  point  on  the  low  pressure  side.  From  the  high  pressure  side 
only  one  characteristic  reaches  the  shock.  The  compatibility  equation 
written  for  that  characteristic,  and  the  three  Rankine-Hugoniot  conditions 
across  the  shock  are  the  necessary  and  sufficient  equations  to  determine 
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the  values  of  P,  u,  and  S  at  the  shock  point  on  the  high  pressure  side, 
plus  the  shock  velocity,  W. 

This  basic  idea  can  be  applied  in  many  different  ways,  with  different 
degrees  of  accuracy  and  sophistication.  Not  all  schemes  fit  equally  well 
into  the  general  logic  of  a  program.  The  present  code  (December  1970) , 
which  was  written  with  a  we 11- separated  treatment  of  every  single  item 
for  a  clearer  understanding,  would  not  easily  accept  the  procedure  sug¬ 
gested  in  Ref.  8  (the  seven  equations,  plus  an  equation  defining  the 
shock  path,  are  solved  simultaneously  by  a  two-level  scheme  which  provides 
an  overall  second  order  accuracy) .  The  procedure  of  the  present  code  is 
conceptually  simpler;  however,  the  accuracy  is  only  of  the  first  order 
and  some  iterations  are  required. 

The  computation  proceeds  as  follows.  First,  the  location  of  the 
shock,  x  ,  is  obtained  at  time  t+At  by  integrating  the  equations 

S 

ds 

_J £  =  w 

dt  W 

with  the  first-order  approximation: 

x (t+At) =x  (t)+W(t) At 
s  s 

Then  the  low  pressure  region  is  identified;  the  two  characteristic  and 
the  particle  path  are  traced  back  from  x  (t+At)  into  the  low-pressure 
region.  Their  locations  at  time  t  are  found  and  values  of  P,  u,  and  S 
are  obtained  by  linear  interpolations  on  the  values  at  mesh  points  at 
time  t.  The  compatibility  equations  are  applied  to  find  the  shock  values 
at  time  t+At  on  the  low  pressure  side.  Finally,  the  Rankine-Hugoniot 
conditions  and  the  compatibility  equation  along  the  characteristic  on  the 
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high  pressure  side  are  solved  by  a  trial-and-error  procedure. 

By  so  doing,  the  shock  computation  depends  only  on  values  at  inter¬ 
ior  points  known  at  time  t.  It  is  thus  explicit,  that  is,  uncoupled  from 
any  other  computation  and  it  can  be  performed  in  an  independent  subroutine. 

In  the  same  spirit,  the  contact  discontinuity  can  be  evaluated  as 
follows.  First,  the  location  xc  of  the  discontinuity  at  time  t+At  is 
determined  by  integrating  the  particle  path  equations 


dx 
_ c 

dt 


=  u 


with  the  first-order  approximation: 


x  (t+At)=x  (t)+u(t)At 
c  c 


Two  characteristics  reach  point  x  from  either  side  of  it.  Their  two 

c 

compatibility  equations  are  solved  simultaneously,  considering  that  P  and 
u  are  the  same  on  both  sides  of  the  discontinuity.  The  values  of  S  on 
either  side  are  carried  over  from  t  to  t+At,  unchanged. 

Once  more,  it  can  be  observed  that  a  more  accurate,  two- level  scheme 
could  be  used;  but,  again,  the  present  procedure  makes  the  computation 
explicit  and  it  is  better  suited  for  the  code  in  its  present  form. 

Finally,  for  the  gradient  discontinuities  two  cases  must  be  distin¬ 
guished.  If  the  discontinuity  propagates  an  expansion,  the  path  of  a 
gradient  discontinuity  is  defined  by 


dx 

dt 


=  u  ±  a 


This  equation,  once  the  proper  choice  of  the  sign  is  made,  can  be 
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integrated  by  a  first  order  approximation  as  the  others  for  x  and  x 

s  c 

above.  The  values  of  P,  u,  and  S  at  both  sides  of  the  discontinuity 
must  be  computed  by  using  information  only  from  the  side  from  which  the 
particles  arrive.  In  fact,  of  the  two  characteristics,  one  is  the  dis¬ 
continuity  itself  and  the  other  carries  information  to  x^  from  one  side 
and  carries  information  away  from  x^  on  the  other  side;  similarly,  the 
entropy  maintains  at  x^  the  value  which  it  has  before  the  particle  reaches 
x^.  On  these  grounds,  it  is  easy  to  work  out  a  routine  for  an  explicit 

treatment  of  x  . 

9 

If  the  discontinuity  propagates  a  compression,  a  shock  will  eventually 
form.  Therefore,  it  is  advisable  to  treat  the  discontinuity  as  a  shock 
of  zero  strength. 


IV.  EXTINCTION  OP  DISCONTINUITIES 

In  an  inviscid  flow,  discontinuities  seldom  disappear  spontaneously. 

A  shock  may  be  extinguished,  for  example,  by  letting  it  reach  a  wall 
capable  of  absorbing  its  full  impact  at  the  right  time,  a  condition  hard 
to  satisfy  in  a  numerical  computation.  A  contact  discontinuity  can  never 
be  eliminated,  in  principle. 

Although  the  program  is  capable  of  handling  a  great  number  of  dis¬ 
continuities,  there  is  no  reason  for  carrying  them  along  and  letting  them 
multiply  ad  infinitum,  unless  their  role  is  significant.  Provisions  are 
taken  in  the  program  to  eliminate  shock  waves  whose  pressure  ratio  is 
below  a  certain  tolerance,  contact  discontinuities  whose  temperature  ratio 
is  below  a  certain  tolerance,  and  gradient  discontinuities  where  the 
difference  between  the  values  of  dP/5x  on  either  side  is  below  a  certain 
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tolerance.  Tolerances  can  be  prescribed  by  the  user  according  to  his 
needs. 

In  any  case,  once  a  discontinuity  is  eliminated,  the  two  points 
from  either  side  of  it  are  merged  into  a  single  point  and  all  mesh  points 
are  redistributed  uniformly  across  the  new  region.  Physical  values  at 
the  new  mesh  points  are  linearly  interpolated  from  the  data  available 
immediately  prior  to  the  merging. 


V.  INTERACTION  OF  DISCONTINUITIES 

It  is  well-known  that  the  way  discontinuities  interact  with  each 
other  depends  on  local  properties  of  the  flow  in  the  vicinity  of  the  point 
where  such  discontinuities  meet.  Consequently,  a  general  analysis  of 
interactions  can  be  made  and  has  actually  been  made  in  the  forties.  Its 
results  are  now  available  in  textbooks  at  the  graduate  level  and  I  do  not 
consider  necessary  to  repeat  them  here.  A  computational  code  may  take 
care  of  the  interactions  by; 

1.  Determining  that  an  interaction  should  occur  within  the  step 
to  be  taken,  and 

2.  Performing  a  local  study  of  the  flow  field,  in  which  the  minute 
regions  between  discontinuities  are  assumed  to  be  regions  of  uniform  flow. 

The  first  problem  is  easily  solved  since  the  present  code  keeps 
track  of  the  size  of  the  regions  between  discontinuities  and  of  the 
slopes  of  the  discontinuities. 


FIG.  3.  SHOCK  REFLECTIONS  AND  INTERACTIONS 


Figure  3  shows  the  interactions  involving  impinging  shocks  only, 
considered  by  the  code.  Cases  a  and  b  are  reflections  of  shocks  on  rigid, 
moving  walls.  The  values  at  D  are  assumed  equal  to  the  values  at  C.  The 
slope  AB  of  the  reflected  shock  is,  by  trial  and  error,  obtained  to 
satisfy  the  Rankine-Hugoniot  conditions  across  the  shock  at  D,  under  the 
assumption  that  behind  the  shock  the  velocity  of  the  particles  equals 
the  velocity  of  the  wall.  Note  that,  after  reflection,  the  number  of 
regions  of  continuous  flow  is  still  the  same  as  before  reflection. 

Case  c  represents  the  coalescence  of  two  shocks  of  the  same  nature. 

In  this  case  the  code  simply  assumes  that  the  two  shocks  merge.  Region  2 
between  them  is  eliminated,  and  the  computation  proceeds  by  determining 
the  location  and  slope  of  a  shock  separating  region  1  from  region  3,  as 
outlined  in  Section  III.  Any  expansion  wave  produced  in  the  interaction 
is  sufficiently  weak  to  be  computed  by  the  code  for  continuous  regions 
without  additional  sophistication.  The  contact  discontinuity  which  should 
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emerge  from  the  merger  point  is  considered  negligible. 

In  case  d,  two  shocks  of  opposite  nature  impinge  upon  each  other, 
are  refracted  and  a  contact  discontinuity  is  generated.  Again,  the  values 
at  D  and  F  are  assumed  equal  to  the  values  at  C  and  E,  respectively.  The 
slopes  of  the  refracted  shocks  and  of  the  contact  discontinuity,  as  well 
as  the  (uniform)  values  of  the  flow  field  parameters  in  regions  2  and  3 
behind  the  merger  point  are  obtained  by  trial  and  error;  the  Rankine- 
Hugoniot  conditions  across  both  refracted  shocks  and  the  conditions  of 
equal  P  and  u  in  regions  2  and  3  are  used.  Note  that  in  this  case,  after 
reflection,  the  number  of  regions  of  continuous  flow  is  increased  by  one. 


t 


FIG.  4.  SHOCK-CONTACT  DISCONTINUITY  INTERACTIONS 


A  contact  discontinuity  cannot  reach  a  rigid  wall.  It  can  interact 
with  a  shock  in  two  ways,  represented  in  Fig.  4.  The  different  pattern 
is  a  consequence  of  the  ratio  of  densities  across  the  contact  discon¬ 
tinuity  before  it  crosses  the  shock.  Pattern  (a)  (refracted  shock, 
refracted  contact  discontinuity  and  a  reflected  shock)  occurs  if  the 
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density  in  region  1  is  greater  than  the  density  in  region  2.  In  the 
opposite  instance,  pattern  (b)  occurs  (no  reflected  shock) .  The  handling 
of  such  cases  follows  the  outline  for  shock  interaction  very  closely. 


VI.  GENERATION  OF  DISCONTINUITIES 

In  the  present  code,  which  does  not  consider  layers  of  different 
gases,  contact  discontinuities  can  only  be  generated  by  shock  interaction, 
as  outlined  in  Section  IV.  Shocks  are  generated  by  abrupt  changes  in  the 
speed  of  a  piston,  or  by  coalescence  of  characteristics.  In  the  majority 
of  problems,  the  time  and  place  where  coalescence  occurs  are  not  directly 
related  to  typical  geometrical  features  of  the  duct  or,  if  moving  pistons 
exist,  to  their  accelerations.  A  method  to  find  the  birthplace  of  an 
imbedded  shock  by  analyzing  the  shape  of  the  pressure  distribution  has 
been  outlined  in  Ref.  7;  it  is  applied  in  the  present  code  and,  in  all 
numerical  experiments  performed  so  far,  has  proved  to  be  efficient. 

If  the  piston  velocity  changes  abruptly  at  time  t,  producing  a 
sudden  compression,  a  shock  originates  at  the  piston  itself  at  time  t. 

In  this  case,  a  new  region  of  continuous  flow  behind  the  shock  must  be 

generated.  The  procedure  is  easily  established  by  introducing  a  ficti¬ 
tious  shock  of  zero  strength  impinging  on  the  piston  at  time  t  and  con¬ 
sidering  the  new  shock  as  a  reflection  of  the  former  one  (case  a  or  b  of 

Fig.  3). 

Finally,  the  piston  acceleration  may  change  abruptly  at  time  t,  pro¬ 
ducing  a  compression.  In  this  case,  the  flow  velocity  remains  continuous 
until  a  shock  forms  by  coalescence  of  characteristics.  The  derivatives 
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of  velocity  and  pressure,  though,  are  discontinuous  along  the  character¬ 
istic  issuing  from  the  piston  at  time  t;  in  other  words,  the  shock  will 
eventually  form  along  a  gradient  discontinuity.  One  may  use  the  technique 
mentioned  above  to  find  imbedded  shocks  in  this  case.  However,  unless 
the  characteristic  itself  is  considered  as  a  boundary  between  two  regions, 
the  calculation  may  become  affected  by  gross  inaccuracies  long  before  the 
shock  has  a  chance  to  be  detected.  It  is  much  better,  thus,  to  locate 
the  point  of  discontinuous  acceleration  in  the  piston  and  to  consider  a 
shock  issuing  from  that  point,  by  proceeding  as  explained  above  for  the 
case  of  discontinuous  piston  velocity.  Since  the  velocity  is  not  discon¬ 
tinuous,  the  shock  will  automatically  begin  with  zero  strength  and  proceed 
for  awhile  along  a  characteristic.  The  advantage  in  so  doing  is  that  the 
characteristic  will  be  singled  out  from  its  inception  as  a  boundary  be¬ 
tween  two  continuous  regions,  and  the  computation  along  the  characteristic 
will  have  all  the  features  of  a  shock  computation.  When  coalescence 
occurs,  the  yet  infinitely  weak  shock  will  be  able  to  pick  up  strength 
and  the  neighboring  points  will  always  be  computed  without  carrying  finite 
differences  across  discontinuities  of  any  type.  An  example  of  the 
improvement  in  accuracy  obtained  by  fitting  a  shock  at  a  point  of  discon¬ 
tinuous  acceleration  will  be  given  in  Section  IX. 

It  should  be  noted  that  proper  identification  of  the  starting  point 
of  a  shock  is  imperative  to  obtain  accurate  results  in  the  adjoining 
regions  of  continuous  flow.  If  the  shock  is  fitted  too  far  from  the 
place  where  it  would  actually  originate,  it  will  not  fulfill  its  mission; 
another  shock  will  tend  to  form  in  the  right  place,  with  the  consequent 
loss  of  accuracy  exposed  in  Section  I.  An  example  of  this  kind  of 
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inaccuracy  is  given  in  Ref.  9,  whose  Figs.  27  and  28  are  reproduced  here 
in  Fig.  5.  Figure  28  shows  (in  heavy  lines)  a  piston  path  and  an  assumed 
shock  path,  whereas  the  actual  shock  should  be  the  light  line  issuing 
from  x=3.0948,  t=1.5.  Figure  27  shows  successive  pressure  distributions 
in  the  region  between  the  piston  and  the  assumed  shock.  Eventually,  (at 
t=2.23)  the  real  shock  overtakes  the  assumed  shock  and  the  accuracy 
improves.  However,  the  transient  is  poorly  described  (note  all  the  dis¬ 
turbing  oscillations  which  a  well-fitted  shock  would  eliminate) .  Note 
also  that  the  number  of  mesh  points  (not  mentioned  in  Ref.  9)  seems  to  be 
quite  large,  whereas  the  method  of  the  present  paper  handles  the  same 
problem  with  only  three  points  between  piston  and  shock. 

If  the  shock  were  fitted  too  close  to  the  origin  of  compression  (the 
piston) ,  the  shock  could  not  function  properly  as  a  relieving  mechanism 
for  the  pressure  surge,  and  again,  the  pressure  distribution  would  become 
accordion-pleated  until  the  shock  had  a  chance  to  move  into  its  right 
location.  All  these  inaccuracies  become  very  important  and,  quite  often, 
catastrophic  in  problems  involving  multiple  shocks  and  their  interactions. 


VII.  TREATMENT  OF  CONTINUOUS  REGIONS 

Having  thus  described  briefly  how  the  most  relevant  features  of  one¬ 
dimensional  flows  can  be  treated,  it  remains  to  say  how  the  computation 
may  be  performed  in  regions  of  continuous  flow.  I  have  pointed  out  in 
Ref.  2  that  a  second-order  accurate  scheme  for  integrating  Euler's  equa¬ 
tions  is  a  reasonably  good  choice,  and  I  have  expressed  my  belief  that 
the  most  popular  schemes  are  almost  equivalent,  as  far  as  accuracy  is 
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FIG.  5.  COMPUTATION  OF  THE  FLOW  PRODUCED  BY  A  PISTON 
SUDDENLY  SET  INTO  MOTION,  ACCORDING  TO  REF.  9 
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concerned.  The  present  code,  written  in  its  original  form  a  couple  of 
years  ago,  uses  the  scheme  No.  6  of  Ref.  2,  page  42.  Such  a  scheme  has 
the  advantage  of  being  symmetrical;  in  a  problem  where  two  pistons  move 
symmetrically  at  both  ends  of  a  duct,  the  results  inside  are  symmetrical 
as  well.  MacCormack's  scheme  (No.  T of  Ref.  2,  page  42)  would  not  provide 
such  a  symmetry,  except  for  few  of  the  first  significant  figures.  The 
latter,  however,  has  some  advantages  in  coding  simplicity  and,  when 
properly  combined  with  schemes  to  treat  discontinuities,  can  provide 
some  reduction  in  computational  time.  Therefore,  a  new  version  of  the 
code  will  be  provided  eventually,  with  a  faster  working  logic,  exploiting 
the  positive  features  of  MacCormack's  scheme. 

The  equations  of  motion  are  written  for  the  case  of  a  duct  of  vari¬ 
able  cross-section  in  the  form: 

Pt+uPx+Yux+yua  =  0 

u  +uu  +  ?P  =0 

t  X  X 

S, +uS  =  0 
t  x 

where  P  is  the  logarithm  of  pressure,  u  the  velocity,  V  is  the  ratio  of 
pressure  to  density  and  S  is  the  entropy.  All  these  quantities  are  made 
nondimen sional  as  explained  in  Ref.  2,  page  22  or  in  a  similar  manner 
(for  inlet  problems)  by  assuming  Pre£  and  Pre£  as  the  values  of  pressure 
and  density  at  the  entrance  section  of  the  duct.  The  quantity  a  is  the 
logarithmic  derivative  t  f  the  cross-sectional  area  with  respect  to  x.  In 
each  region  of  continuous  flow,  the  equations  are  normalized  by  assuming 
new  variables  x  and  T, 
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x-b 

c-b 


T  =  t 

where  b  and  c  are  the  left  and  right  boundary  of  the  region,  respectively. 
The  finite-difference  scheme  is  then  applied  to  the  equations  in  the 
(X,T)  space. 

VIII.  INITIAL  CONDITIONS,  BOUNDARY  CONDITIONS 
AND  DISTRIBUTION  OF  MESH  POINTS 

Everyone  who  has  ever  written  a  computational  code  knows  how  diffi¬ 
cult  it  is  to  make  it  flexible  and  general.  Regardless  of  the  care  taken 
in  the  planning  stage  to  give  the  code  some  generality,  soon  a  request 
for  application  will  show  features  of  the  geometry,  boundary  conditions 
and/or  initial  conditions  which  had  been  overlooked.  One  of  the  advan¬ 
tages  of  the  present  approach  proceeds  from  its  segmentary  structure. 

The  basic  routines  (shock  detection,  shock  fitting,  contact  discontinuity 
fitting,  gradient  discontinuity  fitting,  reflections  and  interaction  of 
discontinuities,  treatment  of  continuous  regions)  are  independent  from 
one  another  and  can  be  used  with  no  changes,  regardless  of  the  number  of 
discontinuities  and  of  the  nature  of  initial  and  boundary  conditions. 
Nevertheless,  I  disclaim  any  attempt  to  generality  in  the  present  stage 
of  evolution  of  the  code.  Originated  two  years  ago  as  a  computational 
tool  for  designing  a  shock-tube  device  known  as  the  Slingshot^,  the 
code  was  limited  to  a  constant  area  duct,  bounded  by  two  moving  pistons? 
the  code  included  provisions  for  detection  of  imbedded  shock,  shock 
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fitting,  reflection  and  interaction  (without  contact  discontinuities) . 

In  its  present  form,  the  code  can  handle  all  the  features  described  in 
the  preceding  sections.  As  far  as  the  boundaries  are  concerned,  any 
combination  of  the  following  conditions  can  be  handled: 

a)  a  moving  piston, 

b)  an  arbitrary  boundary  (either  fixed  or  moving  with  a  prescribed 
law)  in  a  region  of  invariable  flow, 

c)  an  infinite  capacity  containing  gas  at  rest, 

d)  a  prescribed  supersonic  flow, 

e)  an  infinite  cavity  in  which  the  gas  issues  as  a  jet, 

f)  an  interface  with  an  engine,  whose  characteristics  are  to  be 
prescribed  by  an  additional  program. 

The  initial  conditions,  which  were  limited  to  a  gas  at  rest  in  the 
original  version  of  the  code,  may  now  be: 

a)  gas  at  rest, 

b)  steady  flow  in  a  duct  of  variable  cross-section,  with  or  without 
a  normal  shock,  obtained  automatically  by  prescribing  the  exit  pressure. 

Once  the  boundary  and  initial  conditions  have  been  chosen  by  reading 
in  three  code  integers,  and  supplementing  the  duct  geometry  and  the  law 
of  motion  of  the  external  boundaries,  it  remains  to  choose  the  n imber  of 
mesh  points  in  the  initial  configurations.  The  choice  is  dictated  by  the 
duct  geometry  and  the  initial  distribution  of  parameters.  If  the  gas  is 
initially  at  rest,  the  duct  has  a  constant  area  and  the  limiting  pistons 
move  to  produce  compressions,  the  code  will  immediately  split  the  initial 
region  into  three  parts,  separated  by  shocks,  so  that  the  initial  region 
will  still  contain  gas  at  rest  until  the  two  shocks  interact  with  each 
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other.  In  this  case,  the  number  of  initial  points  may  be  limited  to 
three.  If  the  cross-section  has  a  variable  area,  one  needs  a  certain 
resolution  near  the  throat.  A  minimum  of  30  points  should  be  provided 
and  some  increase  may  be  necessary  if  the  time-dependent  computation 
near  the  throat  experiments  difficulties. 

Once  the  initial  choice  of  mesh  points  has  been  made,  no  further 
intervention  is  necessary.  Every  time  a  new  region  is  created,  it  starts 
with  the  minimum  requirement  of  three  points.  As  the  width  of  the  region 
increases,  the  number  of  mesh  points  is  automatically  increased  to  main¬ 
tain  an  adequate  resolution.  Conversely,  if  the  width  of  a  region  de¬ 
creases,  the  number  of  mesh  points  is  reduced  until  the  minimum  of  three 
is  eventually  reached  (an  attempt  to  further  reduce  the  number  of  mesh 
points  results  in  detecting  a  reflection  or  an  interaction,  as  explained 
in  Section  V) . 


IX.  EXAMPLES  AND  DISCUSSION 

The  following  examples  intend  to  show  the  capabilities  of  the  code 
in  its  present  version  and  also  to  illustrate  some  of  the  points  made  in 
the  preceding  sections. 

A.  Cylindrical  Ducts  With  Moving  Pistons 

Let  us  begin  with  the  case  of  a  cylindrical  duct,  bounded  by  two 
moving  pistons,  with  the  gas  initially  at  rest. 

In  Run  No.  2,  whose  results  are  shown  in  Fig.  6,  the  motion  of  the 
left  piston  is  defined  by: 

=  1  -  cos  nt/2  t  <  2 

%  - 2  * »  2 
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IN  A  CYLINDRICAL  DUCT  IMPROVED  COMPUTATION 


The  motion  of  the  right  piston  is  defined  by: 

x  =  3+h  cos  TTt/2  t  <  2 

c 

x  =  2.5  t  >  2 

c 

A  single  region  is  assumed  initially,  in  which  49  points  are  computed 
(including  the  boundary  points) .  No  provisions  are  made  to  define  ficti¬ 
tious  shocks  issuing  from  the  pistons  at  t=0.  Since  the  acceleration  of 
both  pistons  is  different  from  zero  at  t=0,  we  must  expect  two  gradient 
discontinuities  to  propagate  along  characteristics  issuing  from  x=0,  t=0 
and  x=3.5,  t=0.  Both  discontinuities  evolve  into  shocks.  According  to 
Eq.  (39)  of  Ref.  7,  the  shock  produced  by  the  left  piston  should  begin 
at  t=.4,  x=.4732,  the  shock  produced  by  the  right  piston  should  begin  at 
t=.8,  x=2. 5536.  As  pointed  out  in  the  preceding  section,  some  inaccura¬ 
cies  are  bound  to  appear  before  the  shocks  are  formed;  the  shocks  may  not 
be  detected  at  the  right  places  and  the  errors  generated  in  the  initial 
phase  of  the  computation  will  disturb  the  results  in  general.  These 
effects  are  evident  in  Fig.  6.  The  heavy  solid  lines  are  the  piston  and 
shock  trajectories;  contact  discontinuities  are  represented  by  heavy 
broken  lines;  the  lighter  lines  are  lines  of  constant  velocity.  The  two 
circles  indicate  the  theoretical  origins  of  the  two  initial  shocks.  It 
is  clear  that  neither  one  of  the  shocks  is  detected  properly  and  that 
oscillations  tend  to  appear.  The  shocks  not  being  detected  properly, 
their  path  is  slightly  incorrect  (the  right  one  is  displaced  to  the  right 
by  about  .05) ,  their  intersection  does  not  occur  at  the  right  place  and 
the  following  pattern  is  incorrect.  Nevertheless,  the  pattern  is  quali¬ 
tatively  significant,  showing  shock  interactions  and  reflections  at 
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infinitum.  Note  the  disappearance  of  the  contact  discontinuity  as  soon 
as  the  density  jump  across  it  becomes  insignificant.  The  most  important 
conclusion  to  be  drawn  from  this  computation  is  the  waste  of  49  mesh 
points  to  produce  an  imperfect  result. 

The  same  case  is  recomputed  as  Run  No.  3,  shown  in  Fig.  7.  Now  the 
gradient  discontinuities  are  explicitly  taken  into  account  as  dividers 
between  regions,  so  that  from  the  very  beginning  the  computational  field 
is  split  into  three  regions.  The  initial  number  of  mesh  points  is  now 
only  3.  No  oscillations  appear,  the  coalescence  of  characteristics 
(represented  by  constant  velocity  lines  in  the  original  simple  waves) 
occurs  properly,  and  the  shocks  start  gathering  strength  (and,  conse¬ 
quently,  bending  their  paths)  at  the  theoretically  predicted  points.  For 
this  run,  the  tolerance  on  temperature  ratio  across  a  contact  discon¬ 
tinuity  was  taken  so  small  that  no  contact  discontinuity  is  eliminated. 

In  Run  No.  6,  the  motion  of  the  left  piston  is  defined  by 
xb  =  2t  (t<.5) 

xfa  =  -  . 25+t ( 3— t )  ( . 5<t<l . 5) 

xb  =  2  (t>l. 5) 

and  the  motion  of  the  right  piston  is  defined  by 


X  = 
c 

2.5  +  cos  nt/2 

(t<l) 

X  = 

c 

2. 5+  (tt/2)  (2. 5-4t+l.  5t8) 

(l<t<4/3) 

X  = 

2.23833 

(t>4/3) 

The  gas  is  again  initially  at  rest  in  a  cylindrical  duct.  In  this  case, 
the  left  piston  starts  moving  impulsively.  A  shock,  binding  two  regions 
of  uniform  flow,  originates  at  t=0,  x=0  and  moves  at  a  constant  speed 
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until  it  is  overcome  by  the  first  expansion  due  to  the  braking  of  the 
left  piston.  The  computed  results  are  shown  in  Fig.  8,  where  only  piston 
paths,  shocks  and  contact  discontinuities  are  plotted,  for  greater  clarity. 

B.  Flow  in  Ducts  of  Variable  Cross-section 

The  next  experiment  (Run  No.  5)  is  made  on  a  Laval  nozzle,  connect¬ 
ing  an  infinite  capacity  with  an  ambient  where  the  pressure  is  .7  the 
stagnation  pressure.  The  flow  starts  from  a  state  of  rest,  as  if  a  dia¬ 
phragm  were  suddenly  removed  at  the  exit  section.  The  geometry  of  the 
nozzle  is  defined  by 

A  =  x/2  +  1/x  (0  <  x  <  3) 

The  results  (isobars)  are  plotted  in  Fig.  9.  It  can  be  seen,  from  the 
figure,  how  a  shock  builds  up  in  the  divergent  section  of  the  nozzle  and 
develops  until  it  reaches  a  steady  condition.  The  values  at  t=35  (not 
shown  in  the  graph)  are  accurately  representing  a  steady  state. 

Suppose  now  that  the  steady  state,  corresponding  to  the  same  Laval 
nozzle  and  the  same  pressure  ratio,  is  the  initial  state  and  the  pressure 
in  the  exit  chamber  is  quickly  raised  to  the  stagnation  value.  Obviously, 
the  final  state  should  be  a  state  of  rest  throughout  the  nozzle.  The 
results  of  the  computation  for  this  case  are  shown  in  Fig.  10.  In  Fig. 10, 
(Run  No.  9,  isobars  plot)  it  can  be  seen  how  the  pressure  rise  at  the 
exit  is  propagated  upstream  through  the  subsonic  portion  of  the  flow, 
until  it  pushes  the  shock  into  the  supersonic  region.  At  about  t=9.5, 
the  shock  has  reached  the  vicinity  of  the  mouth  of  the  nozzle,  and  its 
strength  has  decreased  so  much  that  the  shock  is  automatically  eliminated 
by  the  program.  An  expanded  region  is  left  in  the  central  portion  of  the 
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FIG.  9 


FROM  REST  BY  REMOVING  A  DIAPHRAGM  AT  THE  EXIT  SECTION 


nozzle,  but  the  pressure  slowly  rises  to  its  asymptotical  stagnation 
value. 

C .  Propagation  of  Numerical  Errors  and  Test  for  Steadiness 

Fig.  11  is  presented  here  to  demonstrate  how  a  study  of  isobar  plots 
helps  in  finding  programming  mistakes.  The  figure  refers  to  the  same 
case  described  by  Fig.  10.  It  is  clear,  though,  that  at  the  very  begin¬ 
ning  of  the  computation  a  disturbance  originates  at  the  entrance  of  the 
nozzle  and  propagates  downstream.  Note  how  a  numerical  disturbance  prop¬ 
agates  as  a  physical  one.  The  disturbance  seems  to  persist  indefinitely; 
indeed,  if  it  had  existed  at  t=0  only,  the  isobars  in  the  supersonic 
region  should  quickly  come  back  to  their  original,  steady  state  locations 
since  they  are  not  affected  by  the  disturbance  at  the  exit  section  until 
the  shock  reaches  them.  This  is  obviously  not  the  case  in  the  figure; 
all  isobars  settle  down  at  a  location  displaced  to  the  right  of  the  initial 
one.  The  rise  in  pressure  at  the  entrance  mouth  is  also  unrealistic. 

The  error  is  rather  small,  though.  Note  that  the  initial  perturbation 
passes  through  the  shock  but  the  shock  itself  is  practically  unaffected. 

All  these  facts  pinpoint  a  mistake  in  the  boundary  conditions  at  the 
infinite  capacity  (x=0) .  Upon  inspection  of  the  program,  it  was  found 
that  the  Mach  number  of  the  gas  at  rest  had  not  been  defined,  and  a  value 
of  .05,  used  for  other  purposes,  had  been  forgotten  where  the  exact  value 
of  0  should  have  been  stored.  Fig.  10  is  the  result  of  the  computation 
after  correcting  the  mistake. 

In  connection  with  the  problem  of  consistency  of  values  computed 
according  to  the  quasi-one-dimensional  steady  state  with  values  obtained 
by  applying  the  unsteady  code  to  a  set  of  steady  initial  and  boundary 
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FIG.  10.  TRANSITION  FROM  STEADY  FLOW  TO  REST  IN  A  LAVAL  NOZZLE 


conditions,  a  number  of  runs  were  made.  One,  for  example,  consists  of 
starting  from  the  same  initial  conditions  as  in  Run  No.  9,  but  letting 
the  exit  pressure  unchanged.  The  pertinent  isobar  figure  is  not  shown. 

It  consists  of  a  set  of  straight  lines,  parallel  to  the  t-axis,  which 
proves  the  perfect  steadiness  of  the  computation,  despite  the  non-uniform 
uniformity  of  values  along  the  nozzle  and  the  variations  in  cross-sectional 
area  (the  run  was  executed  for  several  hundred  time  steps  with  no  changes 
in  the  outputs) . 

A  test  of  this  kind  should  be  conducted,  prior  to  computing  an 
unsteady  case  for  a  variable  area  nozzle,  in  order  to  establish  what  is 
the  optimum  number  of  mesh  points  to  use.  If  the  mesh  points  are  too 
many,  one  wastes  computational  time;  but  if  they  are  too  few,  the  resolu¬ 
tion  may  be  affected.  In  the  latter  case,  initial  steady  results  would 
not  remain  steady  under  steady  boundary  conditions.  This  is  particularly 
important  in  the  vicinity  of  a  throat,  where  changes  in  curvature  of  the 
walls  could  make  the  time-dependent  evaluation  very  delicate. 

D.  Engine  Surge  Simulation  and  Related  Effects 

An  engine  inlet  is  represented  by  a  duct  whose  cross-sectional  area 
is  defined  by 

A  =  1  +  \ ( x-2 )  2  (0  <  x  <  4) 

The  Mach  Number  of  the  flow  at  the  entry  section  is  equal  to  4,  and  the 
ratio,  P°x#  between  exit  pressure  and  entry  pressure  is  30.  Two 

All  runs  described  under  (b)  have  been  made  using  30  intervals  along 
the  Laval  nozzle. 


37 


I 


different  steady  states  are  possible,  one  with  a  shock  in  the  divergent 
portion  of  the  duct  (at  x=3.0948),  which  is  a  stable  configuration,  and 
another  with  a  shock  in  the  convergent  portion  of  the  duct  (at  x=.9052), 
which  is  an  unstable  configuration. 

Starting  from  the  stable  configuration,  the  pressure  at  the  exit 
section  is  let  to  vary  according  to  the  law 

p  =  p°  (1  +  H  sin  t)  (0  < 

ex  ex 

(run  No.  15) . 

At  t=2,  the  exit  pressure  resumes  its  original 
ensuing  flow  is  shown  in  Fig.  12  (isobar  plot)  where 
of  the  duct  affected  by  the  perturbation  is  shown.  The  shock  is 
initially  pushed  toward  the  throat,  but  then  sucked  back  by  the  rapid 
expansion.  The  interesting  feature  of  this  case  is  the  smooth  way  in 
which  the  sudden  jump  in  pressure  at  %-2  is  negotiated  by  the 
computational  technique,  despite  the  presence  of  discontinuities  in  the 
first  derivatives.  It  must  be  noted,  however,  that  20  nodes  had  to  be 
used  behind  the  shock  to  achieve  this  smooth  result.  With  fewer  nodos 
the  transition  produces  an  oscillating  pattern  which  does  not  damp  out . 
Another  interesting  feature  is  a  slight  oscillation,  appearing  at  t=4.5 
in  the  supersonic  region.  This  is  a  purely  numerical  inaccuracy 
produced  by  an  automatic  change  in  the  number  of  nodal  points.  Note 
how  the  numerical  error  is  rapidly  transmitted  into  and  eliminated  by 
the  subsonic  region. 

Several  other  conputations  refer  to  the  case  in  which  the  initial 
steady  state  describes  the  unstable  configuration.  First,  the  exit 
pressure  variations  defined  by 

Pex  =  Pex  <1+J*sin  2t>  (0  <  t  <  TT/2) 

is  used  (Run  No.  8,  Fig.  13,  velocity  plot).  The  initial  pressure  rise 
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coalesces  into  a  shock,  which  quickly  overcomes  the  main  shock,  pushing  it 
out  of  the  duct.  The  following  expansion  does  not  catch  up  with  the  shock. 

In  a  second  problem  (Run  No.  25,  Fig.  14,  velocity  plot)  the  exit 
pressure  law  is 

p  =  p°  (1+.15  sin  3t)  (0<t<TT/3) 

ex  ex 

that  is,  a  weaker  perturbation,  occurring  in  a  shorter  time.  The  com¬ 
pression  waves  do  not  coalesce  into  a  shock  but,  again,  the  expansion  is 

* 

unable  to  recall  the  main  shock  into  the  duct. 

E.  Instability  of  a  Shock  in  the  Convergent  Section  of  a  Supersonic  Inlet 

The  last  figures  present  the  numerical  description  of  the  response  of 

a  shock  to  a  perturbation  in  the  impinging  supersonic  flow  in  the  inlet 
** 

defined  in  D.  The  free  stream  pressure  is  defined  by 

P  =  ±  .15  sin  3t  (0<t<n/3) 

P  =  0  (t>TT/3) 

j.  15  describes  the  flow  corresponding  to  an  increase  in  pressure. 
From  the  figure  we  see  that  the  initial  compression  is  propagated  through 
the  supersonic  portion  of  the  duct;  the  shock  is  pushed  downstream,  and 
the  compression  moves  into  the  subsonic  flow  until  it  reaches  the  infinite 
capacity  at  the  right  end.  There,  it  is  reflected  back  as  an  expansion. 
This  evolution,  occurring  in  a  time  equal  to  about  tt/6,  is  followed  by  an 
expansion  which,  at  the  right  end  of  the  duct,  reflects  back  as  a  com¬ 
pression.  The  latter  acquires  strength  during  its  upstream  motion  until 

•k 

Both  runs  No.  8  and  No.  25  have  been  made  with  a  total  of  30  points 
over  the  entire  region. 

*  * 

These  runs  have  been  made  with  30  mesh  points  over  the  entire  duct. 
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it  coalesces  into  a  shock.  However,  tho  main  shock  cannot  be  switched 
back  to  its  original  position.  Indeed,  the  upstream  supersonic  flow  is 
quickly  reestablished  into  a  practically  steady  state.  The  shock,  thus, 
occurs  at  a  lower  Mach  number  (made  even  lower  by  the  motion  of  the  shock 
in  the  direction  of  the  flow)  and  keeps  losing  strength  until  it  reaches 
the  throat.  From  then  on,  the  shock  gathers  strength  again  and  eventually 
it  locks  in  the  stable  configuration  mentioned  in  D.  Between  t=3.5  and 
t=6.5,  we  note  a  certain  deterioration  in  the  computation  behind  the  main 
shock.  This  is  only  due  to  lack  of  resolution;  the  number  of  mesh  points 
is  indeed  too  small  (about  a  dozen)  behind  the  secondary  shock.  However, 
strong  compressions  are  absorbed  by  several  extra  shocks,  and  the 
pattern  is  quickly  restored  to  normality.  The  figure  give.?,  thus,  a  good 
example  of  the  usage  of  shocks  to  represent  strong  compressions  without 
creating  instabilities  (a  purely  numerical  device  which  is  much  safer  and 
economical,  and  less  dissipative,  than  artificial  viscosity). 

Fig.  16  (Run  No.  42,  isobars  plot)  describes  the  opposite  case  of  a 
decrease  in  pressure  in  the  free  stream.  The  initial  expansion  makes  the 
shock  move  upstream.  When  the  supersonic  flow  recovers  its  steadiness, 
the  shock  occurs  at  a  higher  Mach  number,  that  is,  with  greater  strength, 
and  it  cannot  be  recalled  back. 

The  last  two  figures.  Figs.  17  and  18  (Run  Nos.  43  and  44,  isobars 
plots) ,  confirm  the  stability  of  the  shock  in  the  divergent  section  of  the 
duct;  the  first  under  a  compressi  followed  by  an  expansion,  the  second 
in  the  opposite  case.  Compare  the  two  figures  to  note  how  faster  the 
compression  propagates  with  respect  to  the  expansion.  Compare  Fig.  15 
and  Fig.  17,  Fig.  16  and  Fig.  18  to  note  how  the  initial  pattern,  which 
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FIG.  17.  RESPONSE  OF  A  STABLE  SHOCK  TO  A  COMPRESSION 
FOLLOWED  BY  AN  EXPANSION  IN  THE  MAIN  STREAM 


0  t  .C  i  x  4 

FIG.  18.  RESPONSE  OF  A  STABLE  SHOCK  TO  AN  EXPANSION 

FOLLOWED  BY  A  COMPRESSION  IN  THE  MAIN  STREAM 
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describes  a  compression  (expansion)  in  the  throat  region  in  a  subsonic 
flow,  is  the  opposite  of  the  pattern  which  describes  a  similar  evolution 
in  a  supersonic  flow. 
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NOTE: 

The  present  paper  is  an  attempt  to  illustrate  the  philosophy  and 
the  basic  ideas  of  a  computational  program.  No  important  detail  has 
been  omitted.  However,  a  minute  description  of  all  intricacies  and  sub¬ 
tleties  of  the  code  would  make  the  paper  three  or  four  times  as  longer. 
The  authoi  believe*  that  too  many  details  would  confuse  the  issues  for 
a  reader  interested  in  the  nature  of  the  approach,  and  would  be  required 
only  by  a  restricted  number  of  readers.  The  author  will  be  glad  to  pro¬ 
vide  additional  explanations  to  anyone  who  so  desires. 


